AD-AQyO  0X9  SYSTEMS  SCIENCE  AND  SOFTWARE  LA  JOLLA  CA  F/6  U/3 

NONLINEAR  GROUND  MOTION  FROM  A  MEGATON  NEAR  SURFACE  NUCLEAR  EXP— ETC<U> 
MAR  60  N  RIMER*  E  J  HALDA*  J  T  CHERRY  F19628-79-C-0003 

UNCLASS  IE  IED  SSS-R-6U-4422 _ _ AFGL-TR-80-0167  NL 


ABC  F'LE  ^ AD  A09 0  0 1  9 


1 


AFGL-TR-80-0167 


NONLINEAR  GROUND  MOTION  FROM  A 
MEGATON  NEAR  SURFACE  NUCLEAR  EXPLOSION 


Norton  Rimer 
Eldon  J.  Halda 
J.  Theodore  Cherry 


Systems,  Science  and  Software 

P,0.  Box  1620 

La  Jolla,  California  90238 


Scientific  Report  No,  1 


March  1980 


•  •  -  -  .  •  >  s 


Vv  <\ 


.  V-’ 

<£a 


*3 


j 

Approved  for  public  release;  distribution  unlimited 


AIR  FORCE  GEOPHYSICS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
UNITED  STATES  AIR  FORCE 
HANSCOM  AFB,  MASSACHUSETTS  01731 


P.  O.  BOX  1620,  LA  JOLLA,  CALIFORNIA  92038,  TELEPHONE  (714)  463-0060 


'80  1 


;  |S 


4 


f\  o 


Qualified  requestors  may  obtain  additional  copies  from  the  Defense 
Documentation  Center.  All  others  should  apply  to  the  National 
Technical  Information  Service. 


.X  ■! 


Unclassified 


SECURITY  CLASSIFICATION  OF  This  PACE  rWbtn  0«a  Entered) 


Q32- 

QPTTJ  NUM4 


REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


J.  R^PORTI  NUM*«3 


AFGI^TR-80^0167 


7 


,J.  00 VT  ACCESSION  NO 


S.  RECIPIENT'S  CATALOG  NUMBER 

gb-fMO  QL1  - - 


& 


4.  title  'end  Subtitle) 


NONLINEAR  GROUND  MOTION  PROM  A  MEGATON 
NEAR  SURFACE  NUCLEAR  EXPLOSION* 

*■  * .  .  - . 7 


S.  TYPE  OF  report's 

Scientific 


PERIOD  COVERED 

—  1 


1.  AU-TmOR/aJ  j . —  ■  -  -  W  /  .  — 

;  Norton^imer  /  Eldoi^/Haldal  Ted  Cherry  Q*} 

IfSiZATlON  NAME  ANO  AOORESS 


PtRFOAMiMO  QWS.  REPQR-  f  mBER 

Ml^SS=B=M=AA2Z  _ 


Fl9628-79-C-0J0t^ 


Systems,  Science  and  Software 
P.O.  Box  1620 
La  Jolla,  CA  92039 


10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  A  WORK  UNIT  NUMBERS 


6110 2 F 
tT3TffiS2  AD 


II.  CONTROLLING  office  NAME  ANO  AOORESS 

Air  Force  Geophysics  Laboratory 
Hanscom  AFB,  Massachusetts  01731 
Monitor/K.  C.  Thomson/LWH 


TrrrunT  naTT 


Mar*b  L980 

M  HC  R  ’J  f  MU 

100 


-fWJRSB 


z 


u.  MONITORING  AGCnCy  name  A  AODAESSi'I/  dltloront  from  Controlling  Ottieo) 

r/C 


IS.  SECURITY  CLASS,  (at  thtt  report) 

Unclassified 


Zf-ir/ 


IS*,  declassification.  DOWNGRADING 

schedule 


is.  DISTRIBUTION  STATEMENT  Co/  title  Report) 


Approved  for  public  release;  distribution  unlimited. 


17.  Distribution  statement  (at  (A*  ebetrect  ant  nod  tn  Bloc*  10.  it  dltlerent  train  R*pori; 


■  t.  supplementary  notes 


It.  KEY  WQROS  (Continue  an  reretee  tide  it  nocot torr  and  Homily  by  bloc *  number) 

Ground  coupling;  Airblast;  Shock  wave;  Ground  motion;  Crater; 
Finite  difference  calculations;  Shear  failure;  Tension  failure. 


JO.  ABSTRACT  (Continue  an  reeette  tide  II  neceeeory  end  Identity  by  bloc *  number) 

Presented  here  are  finite  difference  calculations  of  the  ground 
motions  from  a  one-megaton  near-surface  nuclear  explosion  over  a 
saturated  soil.  These  calculations  are  unique  in  that  they  were 
carried  out  to  a  time  of  3.26  sec  so  that  the  motions  could  be 
computed  at  ranges  at  which  the  soil  response  was  linearly  elastic. 
Having  a  detailed  specification  of  the  explosion-induced  ground 
motion  on  cylindrical  monitoring  surfaces  in  the  elastic  regime,  - 


DO 


*  O  Am 

JAN  *1 


1473/1 


EDI  TON  OF  I  NOV  Sf  IS  OBSOLETE 

X  ^  <■ 


Unclassified 


77 


■  •**v  r*  4Ati  Pi/-  a  *9  YwtA  aaGC  A*  •*  0#f«  Entorotf* 


security  classification  of  this  faoe<tfh«i  g«. 


analytical  techniques  from  elasticity  will  be  applied  to 
continue  the  ground  motion  through  layered  earth  models 
out  to  ranges  of  10  to  500  ki-lometersV-y  k > 


The  finite  difference  calculations  were  initialized  at 
0.8  msec  using  feh#  results  of  the  Source  3/5  ground 
coupling  calculation  which  incorporates  a  detailed  model 
of  the  nuclear  device  mass,  materials  and  dimensions. 

The  ground  motions  were  computed  out  to  5.0  msec  in  an 
Eulerian  hydrodynamic  code  and  then  transferred  to  the 
Lagrangian  CRAM  stress  wave  code  containing  appropriate 
constitutive  models,  shear  failure,  tension  failure,  etc. 
A  bowl-shaped  crater  was  calculated  with  a  cratering 
efficiency  of  37*eubi<j— feet  per  ton>of  explosive  energy, 
a  factor  of  5  orj6  smaller  th^n  given  by  empirical 
scaling  formulae/based  on  'thePacific  nuclear  tests. 


-TO 


iieufliTv  classification  of  rms  enter**} 


TABLE  OF  CONTENTS 


Section 


Page 


INTRODUCTION  AND  SUMMARY. 


II.  CONSTITUTIVE  MODELS  AND  MATERIAL  PROPERTIES  .  .  5 

2.1  EQUATION-OF- STATE  .  5 

2.2  SHEAR  FAILURE . 11 

2.3  TENSION  FAILURE . 14 

III.  THE  SOURCE  3/5  GROUND  COUPLING  CALCULATION.  .  .  17 

3.1  INITIAL  CONDITIONS,  RESULTS  AT  800  ySEC  .  .  17 

3.2  SOURCE  3/5  CONVERSION  OF  EQUATION-OF-STATE.  23 

3.3  RESULTS  OF  THE  EULERIAN  CALCULATION  TO 

5  MSEC.  .  .  . . 32 

IV.  NONLINEAR  GROUND  MOTION  OUT  TO  THE  ELASTIC 

RADIUS . 48 

4.1  NUMERICAL  PROCEDURES . 48 

4.2  GROUND  MOTION  RESULTS  .  65 

4.3  CRATER  FORMATION . 8  3 


REFERENCES. 


APPENDIX  A  NONLINEAR  GROUND  MOTION . 97 

APPENDIX  B  ELASTIC  GROUND  MOTIONS . 14  3 


i:  1 : :  t 


LIST  OP  ILLUSTRATIONS 


Figure 

2.1 

2.2 

3.1 

3.2 

3.3 

3.4 

3.5 

3.6 

3.7 

3.8 

3.9 

3.10 

3.11 

3.12 

3.13 

3.14 


Page 


CHEST  equation-of-state  for  saturated  tuff 
showing  constant  energy  lines  (zero  energy 
in  the  plot  corresponds  to  an  ambient  ground 


energy  of  1.342  x  109  ergs/gm) . 7 

Calculated  adiabatic  releases  from  CHEST 
Hugoniot . 8 

Source  3/5  phenomenology  at  800  usee . 18 


Source  3/5  pressure  contours  at  800  usee.  ...  20 

Source  3/5  velocity  vectors  at  800  usee  ....  21 

Tracer  particles  marking  boundary  at  800  usee 
between  air  and  original  ground  material.  ...  22 


Energy  (internal  plus  kinetic)  in  downward- 
moving  ground  material  below  original 
ground  surface  before  and  after  equation-of- 
state  change . 25 

Peak  pressure  versus  range  from  spherically 
symmetric  100  KT  test  calculation . 28 

Pressures  profiles  for  spherically  symmetric 
test  calculations . 29 

Density  profiles  for  spherically  symmetric 

test  calculations . 30 

Velocity  profiles  for  spherically  symmetric 
test  calculations . 31 


Tracer  particles  showing  extent  of  blowoff 
of  original  ground  material  at  1.87  msec.  ...  33 

Tracer  particles  showing  extent  of  blowoff 
of  original  ground  material  at  5.00  msec  ....  34 

Internal  energy  contours  in  ground  material 
at  1.87  msec . . 

Mass  density  contours  in  ground  material 

at  1.87  msec . 37 

Velocity  vectors  at  1.87  msec . 38 


iv 


LIST  OF  ILLUSTRATIONS  (continued) 

Figure  Page 

3.15  Pressure  contours  at  1.87  msec . 39 

3.16  Pressure  contours  at  2.84  msec . 40 

3.17  Pressure  contours  at  3.94  msec . 41 

3.18  Pressure  contours  at  5.00  msec  showing 

main  shock.  . . 42 

3.19  Pressure  contours  at  5.00  msec  out  to  a 

range  of  60  meters . 4  3 

3.20  Pressure  contours  in  ground  at  ranges 

greater  than  60  meters  at  5.00  msec . 44 

3.21  Velocity  vectors  at  5.00  msec . 46 

3.22  Density  contours  at  5.00  msec . 47 

4.1  CRAM  grid  at  10.15  msec . 51 

4.2  Velocity  vectors  at  10.15  msec . 5  2 

4.3  Pressure  contours  at  10.15  msec . 5  3 

4.4  Grid  and  crack  angles  at  26.5  msec . 58 

4.5  Grid  and  crack  angles  at  51.0  msec . 59 

4.6  Locations  of  monitoring  surfaces  in  CRAM 

and  SAGE  (not  to  scale) . 61 

4.7  Complete  CRAM  grid  and  crack  angles  at 

321.91  msec . 62 

4.8  CRAM  grid  and  crack  angles  at  647  msec . 64 

4.9  Pressure  versus  time  at  a  horizontal  range 

of  80  meters  and  a  depth  of  10  meters . 66 

4.10  Horizontal  radial  (R)  and  vertical  components 
(Z)  of  velocity  at  a  horizontal  range  of 

80  meters  and  a  depth  of  10  meters . 67 


4.11  Displacement  components  at  a  horizontal  range 

of  80  meters  and  a  depth  of  10  meters . 68 


v 


LIST  OF  ILLUSTRATIONS  (continued) 


Figure 

Page 

4.12 

Pressure  versus  time  at  range  of  80  meters 
and  depth  of  35  meters  . 

70 

4.13 

Velocities  versus  time  at  range  of  80  meters 
and  depth  of  35  meters  . 

71 

4.14 

Displacements  versus  time  at  range  of  80  meters 
and  depth  of  35  meters . 72 

4.15 

Pressure  versus  time  at  a  range  of  80  meters 
and  a  depth  of  75  meters  . 

73 

4.16 

Velocities  versus  time  at  a  range  of  80 
meters  and  a  depth  of  7  5  meters . 

74 

4.17 

Displacements  versus  time  at  a  range  of  80 
meters  and  a  depth  of  75  meters . 

75 

4.18 

Displacements,  velocities  and  stresses 
versus  time  8  meters  from  the  free  surface 
at  a  range  of  712  meters  . 

77 

4.19 

Displacements,  velocities,  and  stresses 
versus  time  8  meters  from  the  free  surface  at 
range  of  815  meters . 

a 

78 

4.20 

Peak  pressure  in  ground  versus  depth 
directly  below  ground  zero  . 

80 

4.21 

Contours  of  maximum  pressure  penetration 
showing  immediate  source  region . 

8] 

4.22 

Contours  of  maximum  pressure  penetration 
showing  stress  levels  less  than  3  kbars.  .  .  . 

82 

4.23 

Pressure  contours  at  26.52  msec . 

84 

4.24 

Pressure  contours  at  99.39  msec . 

85 

4.25 

Pressure  contours  at  547  msec . 

86 

4.26 

Velocity  vectors  at  26.52  msec  . 

83 

4.27 

Velocity  vectors  at  51.05  msec  . 

89 

4.28 

Velocity  vectors  at  99.39  msec  . 

9Q 

vi 


1  A 


LIST  OF  ILLUSTRATIONS  (continued) 

Figure  Page 

4.29  Velocity  vectors  of  321.91  msec . 91 

4.30  Velocity  vectors  at  647  msec . 92 

4.31a  CRAM  grid  at  1.17  sec  with  ejecta  removed  ...  93 
4.31b  Histogram  of  ejected  material  .  93 


vii 


LIST  OF  TABLES 


Table  Page 

2.1  Elastic  Constants  .  10 

3.1  Total  Energies  in  Calculation  Grid  at 

800  ysec  Before  and  After  Equation-of-State 
Conversion . 26 


I 


4.1  Summary  of  Lagrangian  Calculation 


55 


ACKNOWLEDGEMENTS 


It  is  a  pleasure  to  acknowledge  the  contributions  of 
Doctors  J.  C.  Baker  and  L.  E.  Bailey  and  Mr.  J.  L.  Waddell  in 
the  early  time  Eulerian  stages  of  the  calculation.  Dr.  K.  D. 
Pyatt,  Jr.,  has  been  most  helpful  in  the  analysis  of  the  crater 
related  ground  motions.  The  calculation  from  0.8  to  5.0  msec 
and  subsequent  analyses  was  supported  in  part  from  funds  pro¬ 
vided  by  the  Defense  Nuclear  Agency. 


i.  int:  eduction  and  summary 

The  definition  of  both  the  near  and  far-field  ground 
motion  environment  is  a  critical  input  to  the  design  of  land- 
based  strategic  weapons  systems  that  are  required  to  survive 
a  first  strike  nuclear  attack.  In  order  to  help  this  definition, 
we  present  in  this  report,  a  detailed  computer  calculation  of 
the  explosion  induced  ground  motions  from  a  1-MT  near  surface 
burst. 

The  ground  motions  have  been  computed  to  a  time  of  3.26 
sec.  out  to  ranges  of  800  meters  where  the  material  response 
is  linearly  elastic.  An  important  result  of  the  work  is  the 
output  of  the  calculation  at  monitoring  surfaces  in  the  elastic 
regime.  Analytical  propagation  techniques  from  theoretical 
seismology  can  be  applied  to  this  data  so  that  the  response 
of  geologic  and  structural  configurations  of  interest  can  be 
computed.  In  this  sense  these  results  will  provide  a  resource 
for  the  community  that  can  be  used  for  future  studies.  A  sub¬ 
sequent  report  will  present  the  propagation  of  this  data 
through  earth  models  appropriate  for  MX  valleys. 

Our  calculation  was  begun  using  the  results  of  the 
Source  3/5  calculation  as  an  energy  source.  Source  3/5  was 
an  extensive  calculation  of  the  ground  coupling  from  a  1-MT 
surface  burst  over  wet  tuff  performed  by  S3  for  the  Defense 
Nuclear  Agency  (DNA) .  The  Source  3/5  calculation  began  wi th 
a  highly  detailed  1-MT  source  at  a  58  cm  height  of  burst.  (The 

source  was  scaled  to  the  correct  yield  from  the  previously 
reported  Source  3  (Allen  and  Knowles,  1971;  Allen  and  Schneyer, 
1973).)  The  results  of  the  Source  3/5  calculation  at  800  usee, 
some  of  which  will  be  summarized  in  Section  III  of  this  report) 
served  as  the  starting  point  for  our  calculation. 

The  Source  3/5  calculation  (out  to  800  usee)  used  a  gas 
equation-of-state  for  the  ground  material,  allowing  it  to  have 


1 


higher  compression  than  would  otherwise  be  expected.  After  a 
few  one-dimensional  test  problems  were  run,  it  was  determined 
that  the  pressures  and  velocities  were  reasonable  compared  to 
that  expected  from  a  solid  equation-of-state .  Since  the  prob¬ 
lem  was  to  be  continued  to  late  time  where  solid  effects  be¬ 
come  more  significant,  it  was  decided  to  switch  the  equation- 
of-state  to  a  better  representation  of  a  ground  material.  This 
was  done  by  retaining  the  velocity  field  and  pressure  and  spe¬ 
cific  internal  energy  distributions.  Only  material  density 
was  adjusted.  Test  calculations  indicated  this  would  have  a 
small  perturbation  on  the  ground  shock. 

The  calculation  had  been  run  to  800  usee  with  radiation, 
at  which  time  it  was  determined  that  radiation  was  no  longer 
playing  a  significant  role.  The  calculation  was  continued 
in  the  same  code  (STREAK,  a  two-dimensional  Eulerian  radiation 
hydrodynamic  code)  in  a  hydrodynamic-only  mode  to  a  time  of 
5  msec  using  the  solid  equation-of-state.  At  this  point  the 
calculation  was  transferred  to  CRAM,  a  two-dimensional  La- 
grangian  elastic-plastic  ground  motion  code  containing  con¬ 
stitutive  models  appropriate  for  the  behavior  of  geologic 
materials  over  a  wide  range  of  stresses,  including  a  linking 
of  an  effective  stress  law  for  porous  and  saturated  materials 
with  detailed  tension  and  shear  failure  models. 

As  described  in  Melzer,  et  al.  (1979)  future  MX  sites 
lie  in  alluvial  valleys  which  are  filled  with  a  variety  of 
geological  materials.  The  major  purpose  of  this  calculation 
was  to  study  in  detail  the  far  field  ground  motion  from  a 
1-MT  surface  burst.  For  this  reason,  rather  than  to  calcu¬ 
late  a  particular  MX  site,  we  chose  a  simple  geological  con¬ 
figuration,  a  homogenous,  fully  saturated  ground  material 
having  an  unconfined  strength  of  18.75  bars.  Although  the  con¬ 
stitutive  model  provided  for  a  maximum  strength  of  100  bars, 
applying  the  effective  stress  concept  to  the  saturated  ground 
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reduced  its  effective  maximum  strength  to  under  20  bars. 
Therefore,  the  ground  material  being  modeled  represents  a 
very  weak  saturated  rock  or  competent  soil. 

At  the  time  of  overlay  (5  msec)  the  airblast  being 
modeled  explicitly  in  the  STREAK  code  was  replaced  by  a  time- 
dependent  pressure  boundary  condition  on  the  free  surface  in 
CRAM.  The  Brode  (Brode,  1968)  airblast  solution  for  a  1-MT 
surface  burst  was  used  as  the  boundary  condition  up  until  a 
time  of  approximately  306  msec.  At  later  times,  the  Needham 
(Needham,  et.  al. ,  1975)  representation,  which  is  slightly 
better  at  small  overpressures,  was  used  as  the  free  surface 
boundary  condition. 

A  grid  from  the  two-dimensional  linear  elastic  small 
deformation  SAGE  (Cherry,  et.  al. ,  1974)  code  was  placed 
around  the  outer  boundaries  of  the  CRAM  grid.  This  allowed 
us  to  follow  the  near-field  ground  motion  to  late  times 
(several  seconds)  more  efficiently  without  nonphysical  re¬ 
flections  propagating  back  from  grid  boundaries.  As  a  further 
cost  saver,  we  have  recently  developed  an  absorbing  boundary 
condition  for  SAGE  boundaries  which  provides  a  momentum  trap 
for  almost  all  of  the  incoming  P-waves  and  a  good  deal  of  the 
S-waves. 

The  finite  difference  calculation  has  been  run  out  to 
3.26  seconds.  Data  (time  histories  of  stresses,  velocities, 
and  displacements)  has  been  saved  along  monitoring  surfaces 
on  which  the  motions  are  linear  elastic  in  both  CRAM  and 
SAGE  in  order  to  provide  redundancy  to  verify  our  analytical 
continuation  techniques.  At  2.70  seconds,  any  information 
reflecting  from  the  absorbing  SAGE  boundaries  which  travel 
at  the  P-wave  speed  of  the  medium  (2000  m/sec)  is  calculated 
to  begin  arriving  at- the  outermost  monitoring  surface  in  SAGE 
(a  cylinder  extending  to  a  depth  of  677  meters  and  a  radial 
distance  of  815  meters.  The  other  two  monitoring  cylinders 


are  located  approximately  100  meters  (in  SAGE)  and  200  meters 
(in  CRAM)  closer  in.  If  necessary,  the  finite  difference 
calculation  could  be  run  without  introducing  much  error  to 
4.1  sec.  at  which  time  information  from  the  absorbing  boundaries 
at  the  S-wave  speed  (632  m/sec)  begins  to  arrive  at  the  outer¬ 
most  monitoring  surface. 

Beginning  at  about  0.75  sec,  estimates  of  the  crater 
dimensions  and  crater  volume  were  made  from  the  calculation 
at  regular  intervals  up  to  approximately  1.25  sec.  Each  grid 
element  having  sufficient  momentum  to  escape  from  the  grid 
was  ejected  ballistically  in  order  to  estimate  its  final  lo¬ 
cation  in  the  free  surface.  The  estimates  of  crater  dimensions 
at  various  time  intervals  differed  only  slightly;  representa¬ 
tive  dimensions  at  1.17  sec  were  a  crater  depth  of  approxi¬ 
mately  90  meters  and  a  crater  radius  of  105  meters.  These 
give  a  crater  volume  of  approximately  37  cubic  feet  per  ton 
of  explosive  energy,  in  line  with  the  most  recent  calculations 
of  crater  size.  As  will  be  discussed  later,  the  calculated 
crater  volumes  are  still  smaller  than  empirical  estimates  of 
generic  crater  volumes  for  wet  soft  rock  (100  ft3/ton)  and  wet 
soil  (200  ft^/ton)  based  on  scaling  from  the  Pacific  nuclear 
tests  in  saturated  coral  (Cooper,  1976) . 

The  remainder  of  this  report  is  organized  in  four  sec¬ 
tions.  In  Section  II,  we  discuss  the  constitutive  models 
and  material  properties  used  in  calculating  the  ground  motion. 

In  Section  III  we  summarize  the  results  of  the  Source  3/5 
calculation  at  800  usee.,  discuss  in  some  detail  the  change¬ 
over  to  the  solid  equation-of-state  and  present  results  ob¬ 
tained  using  the  Euler ian  code  out  to  5  msec,  the  time  of  over¬ 
lay  into  CRAM.  Section  IV  begins  with  a  description  of  the 
nonlinear  CRAM  calculation  including  the  approximations 
(boundary  conditions,  regions,  etc.)  used.  This  section  con¬ 
tinues  with  a  presentation  of  the  ground  motion  results 
and  concludes  with  a  discussion  of  the  crater  formation. 
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II.  CONSTITUTIVE  MODELS  AND  MATERIAL  PROPERTIES 

A  fairly  simple  geological  configuration,  a  homogeneous 
fully  saturated  competent  soil,  was  chosen  for  this  calculation. 
This  saturated,  weak  medium,  when  compared  to  a  dry  soil  was 
expected  to  give  less  attenuation  of  the  direct  coupled  ground 
motion  as  it  propagated  away  from  the  source  as  well  as  less 
attenuation  of  the  airblast  induced  ground  motion.  By  using 
a  completely  saturated  soil,  we  avoid  the  modeling  approxi¬ 
mations  inherent  in  crushing  up  air-filled  voids  in  an  asym¬ 
metric  environment  but  not  at  the  expense  of  any  generality 
in  the  far  field  solutions  which  are  of  prime  interest  here. 

A  homogeneous  media  was  chosen  rather  than  a  layered  media  in 
which  reflections  could  obscure  the  main  features  of  the  pro¬ 
pagation. 

The  constitutive  models  used  to  predict  the  nonlinear 
behavior  of  the  ground  material  may  be  separated  into  three 
basic  parts,  an  equation-of- state  to  describe  the  thermo¬ 
dynamic  state  of  the  material  (pressure  as  a  function  of  energy 
and  density) ,  a  shear  failure  model  to  account  for  the  material 
strength  of  the  rock,  and  a  tension  failure  model  to  account 
for  the  material  strength  of  the  rock,  and  a  tension  failure 
model  to  account  for  the  opening  and  closing  of  tensile  cracks 
in  the  rock.  A  complete  description  of  these  constitutive 
models  may  be  found  in  Cherry,  et.  al. ,  1975. 

2.1  EQUAT I ON-OF- STATE 

CHEST  is  a  chemical  equilibrium  equation-of-state  for 
saturated  tuff  which  was  developed  at  S3  (Laird,  1976)  espe¬ 
cially  for  use  with  hydrodynamic  codes  to  study  nuclear  ex¬ 
plosion  phenomenology.  This  tabular  equation-of-state  gives 
pressure  as  a  function  of  energy  and  specific  volume  and 
accurately  models  the  rock  behavior  in  pressure  regimes  from 
tens  of  megabars  down  to  a  few  tenths  of  a  bar.  CHEST, 


plotted  in  Figure  2.1,  coupled  an  elaborate  chemical  equilibrium 

treatment  with  steam  tables  and  bulk  modulus  data  to  provide  one 

consistent  equation-of-state  over  this  entire  energy-specific 

volume  space.  For  calculational  reasons,  an  ambient  ground 
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energy  of  1.342  x  10  ergs/gm  at  one  atmosphere  and  98  degrees 
Kelvin  was  defined  as  zero  energy  in  the  finite  difference 
code.  The  constant  energy  lines  shown  in  Figure  2.1  are  labeled 
relative  to  this  ambient  energy. 

The  CHEST  equation-of-state  covers  regions  where  the 
tuff  is  condensed  and  where  it  is  vaporized.  The  curves  for 
energy  greater  than  12.7  x  10^  erg/gm  in  Figure  2.1,  for  example, 
correspond  to  tuff  that  is  vaporized  and  must  be  treated  as 
a  gas  of  chemically  reacting  constituents.  For  smaller  energy 
densities,  the  curves  lie  in  a  region  where  the  tuff  is  in 
mixed  phase,  partly  vapor  and  partly  condensed.  The  tuff 
mixed  phase  and  tuff  vapor  regions  are  marked  "I”  in  Figure 
2.1.  In  the  region  marked  "XI,"  tuff  is  solid,  but  the  free 
water  content  is  part  liquid,  part  vapor.  Thus,  the  curve 
bounding  region  II  is  the  steam  dome  for  the  water  content  of 
the  tuff.  The  narrow  region  marked  III  in  the  figure  is  the 
domain  in  which  the  tuff  is  solid  and  water  is  liquid.  The 
nearly  vertical  curves  arising  from  the  left  boundary  of  the 
steam  dome  are  lines  of  constant  energy.  The  experimentally 
observed  Hugoniot  for  wet  tuff  is  nearly  parallel  to  these 
curves  up  to  pressures  of  the  order  of  100  kbar. 

The  CHEST  table  is  designed  to  be  used  with  a  high 
speed  computer  interpolation  routine.  Figure  2.2  shows 
adiabatic  releases  computed  from  the  CHEST  Hugoniot  for  wet 
tuff.  This  model  contained  a  total  water  content  of  23.66 
percent  by  weight  (including  bound  water) .  Because  of  the 
importance  of  a  detailed  treatment  of  the  water-rock  mixed 
phase  and  vapor  regions  in  particular,  we  chose  to  use  the 
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Figure  2.1  CHEST  equation-of-state  for  saturated  tuff 
showing  constant  energy  lines  (zero  energy 
in  the  plot  corresponds  to  an  ambient  ground 
energy  of  1.  342  x  *109  ercs/gm)  . 
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CHEST  equation-of-state  of  Figures  2.1  and  2.2  to  model  our 
fully  saturated  ground  material.  CHEST  tuff  has  a  bulk 
density  and  a  bulk  modulus  comparable  to  saturated  clays  or 
silty  sands  found  below  the  water  table  in  valleys  of  interest 
for  MX  sites  (see  Melzer,  et  al . ,  1979).  Table  2.1  gives 
the  elastic  constants  used  in  our  model. 

Since  the  finite  difference  calculation  must  be  run 
to  very  late  times  where  material  behavior  is  elastic  over 
most  of  the  grid  of  interest,  CHEST  was  replaced  in  the  cal¬ 
culation  at  those  times  and  locations  by  a  polynomial  fitted 
to  the  CHEST  Hugoniot  at  these  lower  pressures. 


TABLE  2.1  ELASTIC  CONSTANTS 


Bulk  density 

s 

1.9445 

gm/cm 

Longitudinal  wave 

speed  = 

2000 

m/sec 

Shear  wave  speed 

= 

632 

m/sec 

Bulk  modulus 

= 

67.25 

kb 

Shear  modulus 

7.76 

kb 

Poisson's  ratio 

= 

0.445 
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2.2 


SHEAR  FAILURE 


Cherry  and  Peterson  (1970)  have  identified  the  depen¬ 
dence  of  material  strength  on  the  third  deviatoric  invariant 
and  shown  the  improvement  in  the  definition  of  material 
strength,  Y,  when  pressure,  P,  is  replaced  by  P.  Y,  P  and 
P  are  obtained  from  the  stress  invariants,  ,  J£,  as 
follows : 


Y  =  (3J£)1/2 


Here,  P  is  the  pressure  component  with  the  overburden 
pressure,  Pg,  added.  J^,  J2  and  are  the  first,  second 
deviatoric  and  third  deviatoric  stress  invariants.  If  0^.1' 
a22'  a33  are  Pri-nciPal  stresses,  then 

°11  +  a22  +  °33 

(gu  +  P)2  +  (q22  *  ^  +  <°33  +  P)2 

2 


J3 


(011  +  P*  (a22  +  ^°33  +  P)  ’ 


Note  that  when  the  intermediate  principal  stress, 
a 22 f  is  equal  to  either  the  maximum,  a^,  or  njinimum, 
principal  stress,  then 


and 


Y  =  |axl  -  a22| 

P  =  -  PH  +  q22) 
2 


11 


The  failure  surface  for  the  ground  material  in  a  dry 
state  is  modified  to  account  for  the  presence  of  water  using 
"effective  stress"  principles.  These  state  that  the  effective 
stress,  peff'  given  by 


wnere  if  the  fluid  or  pore  pressure,  should  be  used  in 
defining  the  dependence  of  material  strength  on  the  stress 
state.  Thus  our  failure  surface  was  chosen  (at  low  temperature 
or  internal  energy)  to  be 

Y  =  15.0  +  0.4  P  (bars) 

ef  f 

with  a  maximum  strength  of  100  bars. 

For  our  fully  saturated  material,  we  assumed  that  the 
pore  pressure  was  identically  equal  to  the  mean  stress  in 
the  soil.  Thus  the  effective  stress  degenerates  to 


so  that  the  effective  maximum  strength  is  less  than  20  bars. 

A  further  modification  was  made  to  the  failure  surface 
to  account  for  loss  of  strength  at  high  temperatures  (internal 
energy) .  The  calculated  strength  Y  was  reduced  by  a  multi¬ 
plicative  factor 


where  e  is  the  specific  internal  energy  and  e^  the  melt  energy, 
chosen  to  be  2.05  x  10^  ergs/gm.  Once  a  material  is  melted 
it  is  assumed  to  have  no  strength  even  upon  condensation. 


Hooke's  law  is 
the  stress  deviators, 


used  to  obtain  an  initial  estimate  of 
i  .  e .  , 


+  2u 


At 
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where  and  are  the  values  of  the  stress  deviator  at 

time  t  +  At  and  t,  respectively. 

y  is  the  shear  modulus 

e^  is  the  strain  rate  deviator 

At  is  the  time  increment. 

Shear  failure  occurs  if  the  material  strength  evalu¬ 
ated  at  P  in  our  calculation)  is  exceeded,  i.e.,  if 

A 

J  >  Y  (P) 

where 


J  = 


3  /  "2  22  -2 

7  [  S11  S22  S33 


£  n-t-1  l(SllS22^33  \ 

P  =  p  _  - 2 - / 


AAA 


1/2 
1/3 


Pn+1  includes  the  overburden  pressure  and  we  have 

A 

omitted  the  n+1  superscripts  in  S. .. 

J 

If  J  >  Y  (P)  ,  then  adjustment  of  the  stress  deviators, 
is  required.  We  assume  that 

_n+l  _  ~n+l 
Sij  "  aSij 

where  is  the  adjusted  value  of  the  stress  deviator  and 


Y  ( 


P)  +  T  ( 


sJJ.'  1/3 


a  = 


1~2^ 
7  '  2 - 1 


AAA 


1/3 


^  Y  _ 

b  =  gp-  (evaluated  at  P )  . 


The  equation  for  a  is  obtained  by  approximating  the 

A 

strength  function  at  p  with  a  first  order  Taylor  series  and 


assuming  that  no  adjustment  in  P  occurs  during  shear 
failure. 


2.3  TENSION  FAILURE 

Tension  failure  is  assumed  to  occur  in  the  element  if 
a  principal  stress  is  greater  than  zero  and  if  shear  failure 
has  ever  taken  place.  We  then  apply  the  tension  failure  model 
proposed  by  Maenchen  and  Sack  (1964)  and  introduce  an  inelastic 
strain  normal  to  the  crack.  This  inelastic  strain  is  just 
sufficient  to  zero  the  tensile  stress. 

For  example,  if  ^22  an<^  c33  are  t^ie  t^iree  princi¬ 

pal  stresses  and  if  3^.  is  greater  than  zero,  then  the  adjusted 
stress  (a^,  0 22'  °33^  are  9-*-ven  by 


/V 

4  \ 

4en 

all  =  all 

(k 

+  3  ") 

2  v 

a  =  a22  “ 
22  ^ 

(k 

-  I  ‘J) 

Aen 

A 

2  \ 

a33  =  a33  " 

(k 

"  H 

Aen 

where 
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k  is  the  bulk  modulus,  y  is  the  shear  modulus  and  all  stresses 
include  the  overburden  pressure. 

If  two  of  the  principal  stresses,  say  and 

are  greater  than  zero,  then  the  stress  adjustment  becomes 


a 

a 
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4 

/V 

(k 
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All  the  inelastic  strain  increments  (Ae^,  Ae 22 <  Ae^)  are 
accumulated  on  each  cycle  giving 


,n+l 

_n 

Aeu 

11 

=  E11 

+ 

n+1 

_n 

Ae22 

22 

E22 

+ 

,n+l 

„n 

Ae3  3 

33 

=  E33 

4- 

These  equations  give  the  basic  stress-strain  adjust¬ 
ment  during  tension  failure.  However,  they  apply  equally 
well  for  crack  closure.  If  at  least  one  of  the  total  strains 
(say  E^)  is  greater  than  zero,  then  the  crack  will  open  or 
close  depending  on  the  sign  of  a^.  If  >  0,  then 


Ae 


11 


> 


0 


„n+l  pn 
E11  bll 


and  the  crack  width  increases.  The  inequalities  are  reversed 
if  0^^  <  0  and  the  crack  width  decreases.  Closure  will  con¬ 
tinue  until 


+ 


Ae 


11 


<  0. 


Then 


Ae 


11 

,n+l 

'll 


and  the  crack  is  completely  closed.  When  this  state  is 
achieved,  the  element  is  able  to  support  a  compressive  stress 
in  the  (1,0,0)  direction,  but  is  not  assumed  to  heal. 
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III.  THE  SOURCE  3/5  GROUND  COUPLING  CALCULATION 


In  this  section,  we  briefly  describe  the  early  time 
Source  3/5  calculation  and  detail  the  procedures  and  approx¬ 
imations  made  to  continue  this  calculation  in  an  Eulerian 
code  from  800  usee  (the  initial  time  of  our  ground  motion 
calculation)  out  to  a  time  of  5  msec  when  the  calculation 
was  transferred  to  a  Lagrangian  finite  difference  code. 

3.1  INITIAL  CONDITIONS,  RESULTS  AT  800  ySEC 

Source  3/5,  a  calculation  of  the  energy  source  and 
ground  coupling  out  to  800  usee  from  a  1-MT  near  surface  burst 
(58  cm  height  of  burst) ,  was  scaled  from  the  Pource  3  calcu¬ 
lation  (Allen  and  Knowles,  1971)  in  which  appropriate  masses, 
materials,  and  dimensions  were  used  in  modeling  the  device. 

The  code  used  in  this  calculation  was  STREAK,  a  two-dimensional 
Eulerian  radiation  hydrodynamic  code  used  for  numerous  ground 
coupling  calculations  in  the  past.  In  this  version  of  the 
2D-Vera  family  of  codes,  the  radiation  treatment  was  grey 
(i.e,  no  multifrequency  effects),  nonequilibrium,  time-depen- 
dent  diffusion  with  flux  limiters.  The  hydrodynamic  treat¬ 
ment  was  Eulerian  with  moving  material  boundaries  and  separate 
material  properties  (density,  specific  internal  energy,  and 
velocity)  within  a  mixed  cell.  STREAK  contains  the  ability 
to  package  portions  of  the  space  in  the  calculation  with 
vastly  different  grids  (or  zoning)  and  to  allow  each  package 
to  run  at  its  appropriate  time  step..  The  calculation  was  run 
to  800  usee  (0.800  msec)  with  radiation  on,  at  which  time 
it  was  determined  that  radiation  was  no  longer  playing  a  sig¬ 
nificant  role. 

Figure  3.1  sketches  some  of  the  important  features  of 
the  Source  3/5  solutions  at  800  usee.  First  note  the  presence 
at  zero  time  of  a  cylindrically  shaped,  air-filled  "room" 
located  directly  beneath  the  working  point.  This  room,  having 
38-cm  thick  concrete  walls,  was  intended  to  represent  part  of 
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ure  3.1.  Source  3/5  phenomenology  at  800  usee. 


an  MX  trench  suffering  a  direct  hit  from  a  1-MT  surface 
burst.  At  800  usee,  this  air-filled  room  has  long  since 
been  compressed  by  the  propagating  shock  wave  and  was  no 
longer  being  considered  in  the  Source  3/5  calculation.  How¬ 
ever,  the  initial  presence  of  the  room  has  had  considerable 
influence  on  the  solutions,  particularly  on  the  direct  in¬ 
duced  ground  motion.  The  expected  spherical  shape  of  the 
directly  coupled  shock  wave  has  been  altered  somewhat.  As 
shown  in  the  equal  pressure  contours  of  Figure  3.1,  the  peak 
pressure  is  not  constant  along  this  shock  front.  Figure  3.2 
gives  these  contours  at  800  usee  in  greater  detail  while 
Figure  3.3  shows  the  particle  velocity  vectors  at  that  time 
indicating  the  disturbance  to  the  downward  directed  spherical 
shock.  We  shall  see  that  at  later  times,  the  influence  of 
this  "room"  on  the  shock  is  small. 

Both  Figures  3.2  and  3.3  show  solutions  confined  to 
Package  1  of  the  STREAK  calculation.  Three  separate  packages 
were  used  in  the  calculation;  Package  1,  including  all  of  the 
ground  material  extending  radially  approximately  12  m  (the 
entire  direct  coupled  shock);  Package  2,  monitoring  the  air 
and  fireball  materials;  and  Package  3,  extending  4  meters 
into  the  ground  and  containing  the  ground-fireball  interface 
from  12  out  to  approximately  80  m.  Package  3  was  very  finely 
zoned  vertically  (and  thus  sub-cycled  relative  to  Package  1) 
in  order  to  accurately  model  the  planar  downward  propagating 
shock  generated  by  radiation  deposition  in  the  ground.  This 
planar  shock,  part  of  which  is  sketched  in  Figure  3.1,  has  a 
magnitude  of  approximately  30  kbar  at  800  usee.  Not  shown  is 
the  airblast  front  extending  out  to  about  80  m  with  a  shock 
wave  peak  pressure  of  14  kbar,  and  with  approximately  6  kbar 
pressure  behind  the  shock.  (Near  the  working  point,  fireball 
pressures  are  as  high  as  15  kbar.) 

Figure  3.4  shows  the  locations  at  800  usee  of  the 
tracer  particles  which  mark  the  boundary  at  zero  time  between 


19 


mu  s.ipoxip‘0? 
max  s.tooxio*11 

&  7.00PXIP'01 

SYMBOL  VAU.E 
A  J.MnXto'0'' 

t  09 

B  d.Jpi'vio 

C  r.rrv  kip‘0,> 

0  I.SIMXIO*’0 
E  2.E11XIA*10 
P  3.  PFIXIi.)" 10 
G  x.  jp<  vui‘10 

>i  J>.^?rxio*'° 

I  1.5XHX10*’1 
J  2.5HX10MI 
K.  3.M»IX1o‘" 

L  x.Jorxio’11 


H 


Figure  3.2 


Source  3/5  pressure  contours  at  800  usee. 
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Figure  3.4 


Tracer  particles  marking  boundary  at  800  usee 
between  air  and  original  ground  material. 
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the  ground  material  and  the  air  and  fireball  materials.  At 
800  usee,  hot  ground  materials  directly  above  the  source 
having  very  low  densities  have  been  blown  over  50  meters  out 
of  the  ground.  Pressures  in  this  material  have  been  sketched 
in  Figure  3.1.  The  sonic  line  at  this  time  is  approximately 
at  the  original  ground  surface.  The  velocity  reversal  line, 
defining  the  boundary  between  upward  and  downward  moving 
material,  is  as  shown  in  Figure  3.1  and  is  moving  downward 
with  time. 

3.2  SOURCE  3/5  CONVERSION  OF  EQUATION-OF- STATE 

The  Source  3/5  calculation  was  run  out  to  800  usee 
using  the  EIONX  equation-of-state  computer  subroutines 
(Pyatt,  1966).  EIONX  assumes  local  thermodynamic  equilibrium 
(LTE)  involving  neutral  atoms,  ions,  and  electrons  in 
solving  an  extensive  system  of  coupled  nonlinear  equations 
for  the  concentrations  of  each  constituent.  The  assumptions 
made  in  solving  these  equations,  as  well  as  some  more  basic 
assumption  in  this  gas  equation-of-state,  have  limited 
validity  in  the  shock  waves  calculated  at  800  usee.  The 
EIONX  equation-of-state,  because  of  its  lack  of  solid-like 
models,  overcompressed  material  at  the  shock  front,  but 
modeled  accurately  the  hot  expanded  blowoff  materials  behind 
the  shock. 

A  decision  was  made  to  convert  to  the  CHEST  wet  tuff 
equation-of-state,  described  in  Section  2.1  before  continuing 
the  Source  3/5  calculation  past  800  usee.  A  series  of 
simplified  one-dimensional  spherically  symmetric  blast  wave 
test  calculations  were  performed  in  order  to  evaluate  dif¬ 
ferent  proposals  for  accomplishing  the  equation-of-state 
conversion.  The  final  conversion  accepted  values  of  pressure, 
specific  internal  energy,  and  velocity  from  the  calculations 
made  with  EIONX  and  iterated  on  density  to  get  the  same  pres¬ 
sure  values  with  CHEST.  This  resulted  in  much  less  energy 
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being  present  in  the  shock  front  since  CHEST  reduced  compressed 
densities  significantly.  In  the  hot  expanded  region,  densities 
were  increased.  Figure  3.5  compares  the  distribution  of  energy 
in  downward-moving  ground  material  (upward  directed  motions 
are  not  important  to  the  calculation  at  this  time)  below  the 
original  ground  surface  before  and  after  this  conversion. 
Plotted  are  values  of  internal  energy  (including  radiation 
energy)  plus  kinetic  energy  of  material  moving  downward  per 
unit  of  area  in  cylindrical  strips  of  ground  at  the  radial 
distances  shown. 

Table  3.1  gives  the  integrated  total  energies  before 
and  after  the  conversion.  Over  the  entire  grid,  the  most 
noticeable  change  was  a  large  decrease  in  kinetic  energy  at 
the  shock  front  (approximately  35  percent)  which  resulted 
in  a  decrease  in  total  energy  of  approximately  5  percent. 

Behind  the  shock  (inside  of  6  meters  radially)  internal 
energy  increased  by  about  32  percent. 

We  do  not  believe  that  the  changes  introduced  into  the 
Source  3/5  results  at  0.8  msec  by  the  conversion  of  equations* 
of-state  described  above  have  seriously  affected  the  solutions 
at  later  times.  To  quantify  this,  two  spherically  symmetric 
blast  wave  calculations  were  made  for  a  device  yield  of  100  KT. 
The  first  calculation  used  the  CHEST  equat ion-of-state  to 
describe  the  ground  material  surrounding  the  device.  The 
second  calculation  used  the  EIONX  equation-of-state  until  the 
peak  of  the  ground  shock  had  reached  a  range  of  1350  cm. 

(Peak  pressure  in  the  shock  at  that  range  was  300  kbars,  com¬ 
parable  to  the  Source  3/5  peak  pressure  at  0.8  msec.)  Then, 
EIONX  was  replaced  by  CHEST,  accepting  values  of  pressure, 
velocity,  and  internal  energy  as  was  done  above,  and  iterating 
to  determine  mass  density  from  the  CHEST  equation-of-state. 

This  calculation  was  continued  using  CHEST  out  to  better  than 
twice  the  range  of  1350  cm. 
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5  Energy  (internal  plus  kinetic)  in  downward- 
moving  ground  material  below  oriqinal  ground 
surface  before  and  after  eouation-of-state 
change . 


TABLE  3.1  TOTAL  ENERGIES  IN  CALCULATION  GRID  AT  800  u SEC 
BEFORE  AND  AFTER  EQUATION-OF-STATE  CONVERSION 


Internal  Energy 
(10^®  ergs) 


Kinetic  Energy 
(10^®  ergs) 


Total 

,20 


ergs) 


Total  Grid 


Before 

6.11 

1.93 

8.04 

After 

6.15 

1.46 

7.61 

Inside  6  Meters 

Before 

2.05 

0.50 

2.55 

After 

2.69 

0.46 

3.15 
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Peak  pressure  versus  range  for  the  two  spherically 
symmetric  test  calculations  are  shown  in  Figure  3.6.  The 
second  calculation,  incorporating  the  equation-of-state  con¬ 
version,  at  a  range  of  1350  cm,  gives  peak  pressures  within 
10  percent  of  the  calculation  made  using  CHEST  from  the  very 
beginning.  Figures  3.7  through  3.9  compare  the  flow  fields 
(pressures,  densities,  and  velocities)  for  the  two  calculations 
at  a  time  when  the  shock  waves  have  traveled  out  to  2700  cm, 
twice  the  range  at  which  the  equation-of-state  conversion  was 
made.  The  agreements  between  pressures  and  densities  is 
outstanding  as  is  the  agreement  between  velocities  near  the 
shock  front.  However,  the  velocities  do  differ  behind  the 
shock  front.  This  is  the  region  in  which  extrapolation 
from  the  spherically  symmetric  calculations  to  the  surface 
burst  is  tenuous  at  best.  On  the  whole,  the  results  of  these 
test  calculations  give  us  a  great  deal  of  confidence  in  the 
credibility  of  the  equation-of-state  conversion  over  most 
of  the  active  grid. 

The  EIONX  equation-of-state  was  intended  to  be  used 
at  high  temperatures  and  pressures.  This  gas  equation-of- 
state  had  a  ground  cold  pressure  of  approximately  3.4  kbar, 
which  caused  the  ground  at  large  radii  (smaller,  too,  of 
course,  but  they  had  less  time  to  work)  to  blow  off.  By 
800  usee  the  ground  surface  had  moved  up  to  approximately 
20  cm  and  was  moving  at  a  velocity  of  800  m/sec.  To  correct 
for  this  artificial  blow  off,  the  ground  was  put  back  to 
normal  density,  cold  motionless  values  beyond  a  radius  of 
76  meters  and  its  surface  was  returned  to  zero.  The  shock 
in  the  air  was  at  a  radius  of  'v  83  meters  at  this  time.  Thus 
the  ground  surface  from  76-84  meters  is  not  impacted  with 
the  correct  airblast  shock  values. 

Also,  in  the  equation-of-state  switch,  there  were  a 
few  zones  at  the  outer  edge  of  what  was  accepted  in  the 
ground  which  gave  slightly  anomalous  values.  Due  to  the 
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Figure  3.9  Velocity  profiles  for  spherically 
symmetric  test  calculations. 


sensitive  density  dependence  of  the  CHEST  equation-of-state , 
a  few  zones  at  radii  of  70-75  meters  and  a  few  centimeters 
depth  had  oressures  of  one  bar  bracketed  between  zones  of 
10  kbar  pressures.  It  was  felt  that  this  would  not  be  a 
significant  perturbation  to  the  calculation  of  the  crater 
since  this  is  in  material  which  should  be  ejected  ultimately. 
During  the  course  of  the  calculation  it  was  also  felt  that 
much  of  this  would  be  washed  out  in  running  or  rezones  and 
would  not  effect  the  far-field  ground  motion. 

3.3  RESULTS  OF  THE  EULERIAN  CALCULATION  TO  5  MSEC 

After  the  conversion  to  the  CHEST  equation-of-state, 
the  calculation  was  run  out  to  5  msec  using  the  Eulerian 
STREAK  code  in  hydrodynamics  mode  with  radiation  turned  off. 
At  5  msec,  the  calculation  was  overlaid  into  the  Lagrangian 
CRAM  code  and  strength  effects  included  for  the  first  time. 
Strength  effects  for  the  first  5  msec  of  the  calculation  are 
negligible  relative  to  the  high  pressures  calculated. 

Figures  3.10  and  3.11  show  the  locations  of  tracer 
particles  at  1.87  and  5.0  msec  defining  the  boundary  of 
original  ground  material  which  has  been  blown  off  at  very  low 
densities.  At  1.87  msec  ground  material  has  reached  a  height 
of  approximately  100  meters  while  at  5.0  msec  it  extends 
about  150  meters  into  the  air.  The  radial  extent  is  approx¬ 
imately  75  meters  at  5.0  msec.  For  comparison,  at  5.0  msec 
the  airblast  has  propagated  out  to  a  radial  range  of  approx¬ 
imately  175  meters. 

Figure  3.12  shows  contours  of  equal  internal  energies 
in  material  still  in  the  ground  at  1.87  msec.  Figure  3.13 
shows  mass  density  contours  for  the  same  material.  The  soil 
directly  under  ground  zero  contains  the  highest  internal 
energies  and  has  been  expanded  to  less  than  one  quarter  of 
the  initial  density.  This  material  is  rapidly  blowing  out  of 
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TRACER  PARTICLES 


Figure  3.10  Tracer  particles  showing  extent  of  blowoff 
of  original  ground  material  at  1.87  usee. 
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Figure  3.11  Tracer  particles  showing  extent  of  blowoff 
of  original  ground  material  at  5.00  msec. 
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Figure  3.12  Internal  energy  contours  in  ground  material 
at  1.87  msec. 
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the  ground  at  particle  velocities  greater  than  one  meter  per 
millisecond.  (See  Figure  3.14  for  velocity  vectors  at  this 
time.)  Figures  3.13  and  3.14  show  the  main  shock  wave  at  a 
spherical  radius  of  roughly  19  meters  from  ground  zero.  The 
intersection  of  this  shock  and  the  planar  fireball-induced 
shock  lies  near  the  free  surface  and  at  a  horizontal  range 
of  approximately  16  meters.  The  planar  shock  extends  out 
to  approximately  100  meters  at  1.87  msec. 

Figures  3.12  through  3.14  show  only  the  results 
which  are  included  in  Package  1  of  the  computational  grid. 

Most  of  the  planar  fireball  and  airblast  induced  shocks  have 
been  calculated  using  Package  3  which  is  very  finely  zoned 
vertically.  All  computational  packages  have  been  rezoned 
several  times  in  the  course  of  these  calculations.  We  shall 
see  that  Package  1  dimensions  are  increased  with  time  to  in¬ 
clude  all  of  the  expanding  spherical  shock  wave. 

Figures  3.15  through  3.18  show  peak  pressure  contours 
from  Package  1  at  four  times  from  1.87  to  5.00  msec.  These 
show  the  propagation  of  the  spherical  shock  front  and  the 
decay  in  peak  pressure  with  time.  By  5  msec,  peak  pressures 
have  decreased  to  below  60  kbars  over  most  of  this  shock 
front  (68.6  kbar  is  the  maximum  in  the  grid).  Figure  3.19 
shows  pressure  contours  at  5.00  msec  including  much  of  the 
planar  airblast-induced  shock  out  to  60  meters.  (Note  the 
change  in  contour  scale  from  Figure  3.18.)  Peak  pressures  in 
the  planar  shock  wave  are  between  5  and  10  kbar  or  approxi¬ 
mately  10  percent  of  the  peak  pressures  in  the  main  shock 
wave.  The  minimum  pressure  calculated  was  0.44  kbar  in  this 
region.  Figure  3.20  shows  pressure  contours  in  the  planar 
shock  at  even  larger  Horizontal  ranges  at  5.00  msec.  Pres¬ 
sures  of  0.5  kbar  extend  out  to  160  meters  (airblast  peak  pres¬ 
sures  of  approximately  1.5  kbar  were  experienced  at  the  free 
surface  at  a  range  of  160  meters).  Note  that  these  pressures 
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Figure  3.13  Mass  density  contours  in  ground  material 
at  1.87  msec 
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Figure  3.18  Pressure  contours  at  5.00  nsec  showing  main 
shock . 
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are  quite  large  when  compared  with  an  assumed  shear  strength 
of  less  than  20  bars  for  the  saturated  soil. 

Velocity  vectors  at  5.00  msec  are  shown  in  Figure  3.21. 
The  largest  of  these  as  before  are  directed  out  of  the  ground 
near  the  burst  point.  Velocities  in  the  airblast-induced 
planar  shock  are  considerably  smaller  than  in  and  behind  the 
main  spherical  shock. 

Figure  3.22  shows  mass  density  contours  at  5.00  msec. 
Again,  the  lowest  densities  are  in  the  high  velocity  blow  off 
region  near  ground  zero.  All  material  with  density  less  than 
1.0  gm/cc  is  very  likely  to  be  blown  off  in  time.  Figure  3.11 
gave  the  boundary  of  original  ground  material  at  5.00  msec. 
This  material  with  very  low  densities  has  not  been  included 
in  Figure  3.22  and  will  be  neglected  in  our  calculation. 
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IV.  NONLINEAR  GROUND  MOTION  OUT  TO  THE  ELASTIC  PADIUS 


The  Source  3/5  ground  coupling  calculation  out  to 
5.00  msec  was  made  using  the  two-dimensional  Euler ian  hydro- 
dynamic  STREAK  code.  At  5  msec,  this  calculation  was  trans¬ 
ferred  to  the  two-dimensional  Lagrangian  stress  wave  code, 
CRAM,  and  material  strength  introduced  into  the  calculation 
for  the  first  time.  In  this  section  we  describe  the  numerical 
procedures  used  and  present  the  results  of  this  calculation. 

4.1  NUMERICAL  PROCEDURES 

Both  Eulerian  and  Lagrangian  computer  codes  have  cer¬ 
tain  advantages  and  disadvantages  for  two-dimensional  ground 
motion  calculations.  A  Lagrangian  code,  because  each  compu¬ 
tational  cell  follows  the  motion  of  a  specified  mass  of  mate¬ 
rial,  gives  a  more  accurate  description  of  the  behavior  of 
that  mass  element,  particularly  when  complex  material  models 
are  needed.  Since  it  was  also  necessary  to  follow  the  ground 
motion  out  to  low  stress  ranges  (a  few  bars),  where  the  dif¬ 
fusion  inherent  in  Eulerian  codes  would  be  likely  to  obscure 
the  physical  solutions,  the  Lagrangian  CRAM  code  was  used  for 
our  calculation. 

Our  greatest  difficulties  in  the  course  of  this  cal¬ 
culation  were  in  overcoming  the  large  zone  distortions  inher¬ 
ent  in  the  Lagrangian  approach  to  a  surface  burst  calculation. 
From  the  beginning,  it  was  obvious  that  we  would  be  unable 
to  adequately  treat  the  hot,  high  velocity,  low  density  blow- 
off  materials  near  ground  zero  which  were  described  earlier 
in  this  report.  Zone  distortions  would  be  just  too  large. 
Therefore  we  decided  to  remove  these  materials  from  consider¬ 
ation  as  they  moved  well  above  ground  level.  Similarly,  we 
chose  not  to  continue  numerically  modeling  the  mixture  of  hot 
gases,  air  and  blowoff  material  above  the  air-ground  surface. 
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We  replaced  these  by  using  as  free  surface  boundary  conditions 
the  airblast  solutions  developed  by  Brode  (1968)  and  Needhan' 
(1975).  The  pressure  as  a  function  of  ranae  along  the  free 
surface  was  specified  at  each  computations.!  time  step  at 
first  using  a  numerical  fit  to  the  Brode  solutions.  Later, 
at  approximately  306  msec,  this  fit  was  replaced  by  the 
Needham  computer  subroutine  which  was  thought  to  be  a  more 
accurate  representation  for  low  pressures. 

The  first  step  in  our  calculation  was  to  overlay  the 
Eulerian  solutions  at  5.00  msec  into  a  Lagrangian  grid  suit¬ 
able  for  continuing  the  simulation.  Relatively  fine  zoning 
was  needed  in  the  vertical  direction  near  the  free  surface 
in  order  to  adequately  treat  the  airblast  induced  motions. 

We  chose  a  vertical  mesh  size  of  0.5  meters  at  the  free  sur¬ 
face  and  very  gradually  increased  mesh  size  with  depth. 
Horizontal  zone  sizes  in  the  cratering  region  were  between 
1  and  2  meters  with  the  smallest  mesh  sizes  ahead  of  and  in 
the  region  of  the  propagating  main  shock.  These  mesh  spacir.cs 
can  be  considerably  coarser  than  those  in  the  Eulerian  c FREAK 
grid  at  5.00  (0.20  to  0.60  meters  both  vertically  and  horizon¬ 
tally)  due  to  the  more  accurate  Lacrangien  definition.  The 
solutions  at  5.00  msec  were  first  rezoned  into  a  STREAK  grid 
identical  to  the  proposed  CRAM  mesh.  Thus  the  overlay  of 
masses,  velocities,  densities,  pressures,  and  specific  interna 
energies  from  STREAK  Packages  1  and  3  to  CRAM  was  accomplished 
very  simply  on  a  one-to-one  zone  basis. 

The  CRAM  zone  boundaries  as  defined  at  5.00  msec  were 
identical  to  the  STREAK  tracers  (see  Figure  3.11)  locatmc  the 
ground-air  interface  at  ranees  from  76  meters  outward.  At 
76  meters,  the  ground  surcace  was  displaced  0.13  m  upward. 

In  closer  than  76  meters,  the  hot,  low  density  blowoff  mate¬ 
rials  above  oriamal  ground  surface  were  not  included  m  the 
calculation.  Figure  3.21  indicated  that  these  materia  is  are 
moving  upward  at  velocities  creator  than  several  kilometers 
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per  second.  CRAM  zones  containing  the  hot  melted  and  gaseous 
materials  still  below  initial  ground  surface  were  flagged  to 
indicate  present  or  earlier  melt.  Our  constitutive  model 
did  not  allow  deviatoric  stresses  in  this  melted  material. 

Elsewhere,  strength  was  included  using  the  model  defined  in 
Section  2.2.  Gravity  was  also  included  in  the  momentum 
equation  for  the  first  time  and  all  inactive  zones  given  a 
hydrostatic  overburden  pressure. 

After  a  careful  check  of  the  overlay,  the  Lagrangian 
calculation  was  begun.  Figure  4.1  is  a  plot  of  a  part  of 
the  CRAM  grid  at  a  time  of  10.15  msec,  100  cycles  into  the 
Lagrangian  calculation.  (The  same  plot  is  used  to  show  ten-  ,  t 

sile  crack  orientation;  no  cracks  were  open  at  this  time). 

Figures  4.2  and  4.3  show  velocity  vectors  and  pressure  con¬ 
tours  over  the  same  portion  of  the  grid.  As  in  earlier  plots 
from  STREAK,  the  positive  Z  direction  is  out  of  the  ground. 

Low  density  soil,  moving  upward  near  ground  zero  at  the  largest 
velocities  shown  by  Figure  4.2,  has  resulted  in  very  large 
distortions  in  the  Lagrangian  grid.  These  distortions  will 
become  worse  with  time.  Elsewhere  the  grid  has  not  distorted 
greatly  from  the  overlay  time.  Some  distortion  can  be  seen 
in  Figure  4.1  marking  motion  at  and  behind  the  main  shock 
wave . 

The  calculation  at  this  time  has  been  zoned  only  to  a 
depth  of  97  meters.  To  the  right,  zoning  is  in  place  out  to 
almost  300  meters  (not  all  shown  in  Figures  4.1  -  4.3)  in 
order  to  accommodate  the  airblast  induced  ground  motion.  At 
10.15  msec,  this  motion  is  out  to  approximately  220  meters 
horizontally  More  zones  will  be  added  before  either  the  air- 
blast  arrival  reaches  the  right  of  the  grid  or  the  direct 
induced  ground  motion  reaches  the  bottom  of  the  grid.  Well 
before  these  occur,  the  continuing  distortion  of  the  grid 
near  ground  zero  will  require  rezoning. 
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Figure  4.2.  Velocity  vectors  at  10.15  msec. 
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Figure  4.3. 


Pressure  contours  at  10.15  msec. 


Our  calculation  was  accomplished  in  several  different 
stages,  successive  benchmarks  in  time  in  general  being  sep¬ 
arated  by  a  rezone  of  the  existing  grid  and  an  overlay  into 
a  newer  computational  grid  covering  a  larger  area  of  activity. 
Major  rezones  were  made  at  times  of  18,  54,  126,  306  and  594 
msec.  Table  4.1  summarizes  the  numerical  procedures  applied 
in  each  stage  of  the  calculation,  listing  the  types  of  rezones 
used  and  the  appropriate  zone  sizes  in  the  cratering  region. 
Note  that  the  minimum  vertical  zone  size  in  that  region  is 
also  the  vertical  zone  size  elsewhere  along  the  free  surface. 

Several  types  of  rezones  are  available  in  CRAM.  The 
most  straightforward  of  these  combines  neighboring  zones 
two  for  one.  All  zones  adjacent  to  a  specified  vertical  or 
horizontal  line  must  be  combined  with  their  neighbors  across 
the  specified  line  to  maintain  the  code  algorithm.  This  is 
equivalent  to  removing  a  horizontal  or  vertical  grid  line 
from  the  mesh,  thus  lowering  the  cost  of  the  calculation  both 
by  increasing  the  time  step  when  the  finer  definition  is  no 
longer  required  and  by  decreasing  the  total  number  of  zones 
to  be  calculated. 

It  is  the  "hand  rezoning"  feature  of  CRAM  that  proved 
most  useful  during  this  calculation.  A  series  of  computational 
subroutines  have  been  developed  which  allow  a  distorted  re¬ 
gion  of  the  mesh  to  be  rezoned  into  a  completely  new  set  of 
regular  grid  lines  or  curves.  The  only  limitations  on  this 
new  set  of  lines  is  that  they  be  continuous  with  the  remainder 
of  the  old  grid.  We  have  used  this  feature  to  remove  low 
density  material  which  is  above  initial  ground  surface  and 
continuing  to  blow  out  of  the  around  with  time,  thus  distortina 
the  Lagrangian  grid.  For  example,  if  we  had  chosen  to  hand 
rezone  the  grid  shown  in  Figure  4.1,  wo  would  have  raker,  -ho 
region  near  ground  zero  out  to  20  or  30  meters  horizontally 
and  down  to  about  25  meters  vertically  and  set  up  a  completely 
new  set  of  grid  lines  each  continuous  with  the  appropriate 
line  at  the  region  boundary.  A  reasonable  choice  w<~>uld  have 
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TABLF  4.1  SUMMARY  OP  LAGRANOIAN  r'nLCPLATION 


Time  Rezone  Summary  New  Zone  Sizes  Other  Features 

(msec)  in  100x100  meter 

cratering  region 


1-2  m  Addition  of  material 

0. 5-2.4  m  strength,  aravit/  and 
overburden  for  the 
first  time.  Brode 
airblasi  bound. try 
condition  replaces 
Package  2  of  the 
STREAK  comp  if  a  t lona 1 
grid . 

18  1.  Hand  re  zone  of  IP,  =  1-3  m 

crater  region  and  AZ  =  1-2.4  m 
removal  of  blowoff 
material  above 
ground  zero. 

2.  Remove  every 
other  horizontal 
grid  line  to  depth 
of  24  m. 

3.  Grid  extended 
to  R  =  480  m. 

Z  =  190  m. 


5  Overlay  from  STREAK  AR  = 

to  CRAM.  Zone  out  A Z  = 
to  R  =  288  m  and 
Z  =  97  m 


54  1.  Same  as  above. 

2.  Remove  every  AR  =  2-3  m 

other  horizontal  AZ  =  1-2.4  m 
from  24-55  meters. 

3 .  Remove  every 
other  vertical  line 
out  to  72  meters. 

4.  Final  CRAM  grid 
dimensions  reached 
R  =  633  m ,  Z  - 

509  m. 


126  1.  Same  as  above  IP  - 

2.  General  2  to  1  AZ  - 
rezone  of  CRAM  grid 
removing  every 
other  horizontal 
and  vertical  line 
with  a  few  excep¬ 
tions  . 


4-5  n  Surround  CRAM  by 

2-4  m  SAGE  grid  to 

R  =  3003  m ,  Z  - 
2  8  5  6  n  with  a  b '  r  b  i  n  a 
boundary  coni i t : ; o o . 
Mon  i  tonne  sf.it  i  ons 
activated  in  cram  nd 
SAGE . 
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TABLE  4.1  SUMMARY  OF  LAG RANG IAN  CALCULATION  (Continued) 


Time  Rezone  Summary  New  Zone  Sizes  Other  Features 

(msec)  in  100x100  meter 

cratering  region 


306 

1.  Same  as  above 

AR  =  4-5  m 

AZ  =  2-4  m 

Replace  Broae  solu¬ 
tion  by  Needham  sub¬ 
routine.  Pressure 
boundary  condition 
applied  directly  to 
all  melted  cells. 

594 

1.  General  2  to  1 
rezone  of  CRAM 
grid . 

AR  =  8-10  m 

AZ  =  4-8  m 

1730 

1.  Removal  of  a 

AR  =  10  m 

Motion  frozen  for 

few  lines  which 
control  time  steo 

A  3  =  8  m 

R  <  140m,  Z  <  140m 
since  not  enough  de¬ 
finition  left  for 

crater  reaion. 


been  to  angle  the  free  surface  line  up  to  Z  =  +5  meters  at 
the  axis  of  symmetry,  throwing  away  the  low  density  material 
above  the  new  free  surface  line. 

In  actuality  we  rezoned  at  IS  msec  when  the  free  sur¬ 
face  was  approximately  30  meters  above  initial  ground.  Figure 
4.4  is  a  plot  of  part  of  the  new  grid  at  26.52  msec,  show  me 
the  larger  vertical  zone  sizes  near  the  free  surface  behind 
the  main  shock  (now  at  a  depth  of  about  90  meters)  as  well 
as  the  reordered  grid  near  ground  zero.  Figure  4.5  plots  the 
grid  and  crack  angles  at  51.05  msec.  (Tensile  cracks  can  be 
seen  at  about  100  meters  from  the  axis  of  symmetry.  These 
will  be  discussed  later.)  The  large  zones  near  ground  zero  in¬ 
dicate  very  low  density,  hot  material  down  to  nearly  40  meters 
and  out  to  about  30  meters.  Adjacent  to  the  free  surface  out 
to  about  70  meters,  zones  which  have  had  extensive  radiation 
deposition  until  0.8  msec  are  at  a  very  low  density.  P.e- 
zoning  of  the  grid  was  accomplished  at  54  msec,  including 
removing  this  low  density  material  as  well  as  some  of  "he 
blowoff  material  near  ground  zero. 

The  calculation  was  run  out  to  126  msec,  at  voich 
time  another  rezone  was  needed.  At  this  time,  the  elastic 
radius  appeared  to  be  expanding  very  slowly  if  at  ail.  There¬ 
fore  a  grid  from  the  two-dimensional  linear  elastic  SATE  mme 
(Cherry  and  Halda,  1974)  was  olaced  around  the  existin'?  CRAM 
grid  out  to  a  radial  range  of  3001  meters  and  a  depth  of 
2856  meters  (R  =  633  m,  2  -509  r  defined  the  CRAM-S/V”: 

boundaries) .  The  SAGE  code  contains  an  absorbing  boundary 
treatment  which  effectively  trams  almost  all  of  the  incident 
P-wave  signals  and  a  creat  deal  of  the  incident  S -waves .  f’h  is 
code  is  often  used  at  S  '  to  provide  chcar,  efficient  bounder 
conditions  for  CP. AP*  (it  is  an  order-  :f  muonitude  less  exten¬ 
sive  to  calculate  ground  tm-ion  in  FACT*  rather  tha  r  CRAM.. 

At  126  msec,  monitor  in-,  stat.i  ns  at  which  velocities, 
displacements  and  stresses  wer  e  to  bo-  s  wed  were  :  l  aced  in 


METERS 


SOURCE  3/5 


both  CRAM  and  SAGE  outside  of  the  elastic  radius.  Figure  4.6 
sketches  the  cylindrical  monitoring  surfaces,  located  to  pro¬ 
vide  redundant  surfaces  with  which  to  check  out  the  analvtica 
continuation  methods.  As  a  failsafe  procedure  to  insure 
linear  elastic  behavior  near  the  monitoring  surfaces  in  CRAM , 
zones  within  25  meters  of  these  surfaces  were  required  to 
behave  elastically.  Unfortunately  our  estimate  of  the  depth 
of  yielding  was  in  error.  CRAM  cells  below  the  mandated 
elastic  surface  at  a  465-meter  depth  would  have  yielded  at 
later  times  if  allowed  to.  A  careful  examination  of  the 
stress  deviators  established  that  preventing  plastic  yielding 
was  equivalent  to  the  soil  below  465  meters  out  to  a  horizon¬ 
tal  range  of  approximately  200  meters  having  an  effective 
strength  of  approximately  32  bars.  Since  realistic  MX 
geologies  would  almost  certainly  be  stronger  at  depth,  no 
error  has  been  introduced  into  the  calculational  results. 
Plastic  flow  did  not  occur  near  the  vertical  monitoring  sur¬ 
face  in  the  calculation. 

The  next  benchmark  in  the  calculation  occurred  at 
306  msec.  At  that  time,  a  hand  rezone  was  aaain  made  in  the 
cratering  region.  The  zones  inside  of  the  melt  radius 
(approximately  60  meters)  were  at  extremely  low  densities 

cr 

(some  at  10  gms/cc) ,  outside  of  the  ranoe  of  the  CHEST 
equation-of-state  tables.  Therefore  we  decided  to  disregard 
the  pressures  given  oy  CHEST  for  those  cel’s  and  to  apply 
the  airblast  pressure  to  the  boundary  o *-ho  ev  dens  i  tv  melt 
This  is  equivalent  to  removing  those  cei  Is  from  the  computa¬ 
tion  for  the  rest  of  the  calculation.  At  this  time,  the 
Brode  airblast  solution  was  replaced  us  a  boundary  next;  o', 
by  the  Needham  airblast  solution,  deemed  to  be  mto  o  iccurato 
at  low  pressures.  Figure  4.7  shows  the  com  1  etc  CRAM  uri.i 
at  321.91  msec,  soon  after  this  rezeno.  The  crater  roc  ion 
at  this  time  is  surrounded  by  cells  contain  ing  tens i !e  cracks 


The  last  major  rezone  was  made  at  594  msec.  Figure 
4.8  shows  the  rezoned  grid  at  647  msec.  Note  that  the  grid 
lines  in  the  melted  region  near  ground  zero  are  completely 
decoupled  from  the  rest  of  the  grid.  Only  at  the  boundary 
of  this  melt  is  the  motion  calculated  by  CRAM.  By  1.73  sec¬ 
onds,  the  crater  calculated  using  ballistic  trajectories  has 
remained  virtually  the  same  for  some  time.  Stresses  and 
velocities  are  very  small  as  well.  Therefore  we  felt  reason 
ably  safe  in  freezing  the  motion  near  the  crater  (R  <  140  m 
and  Z  <  140  m)  and  rezoning  away  the  few  lines  controlling 
the  time  step  so  that  the  ground  motions  at  the  monitoring 
surfaces  could  be  calculated  to  later  times  more  cheaply. 

The  calculation  with  this  final  grid  was  run  out  to  3.26  sec 
The  plots  of  velocities,  displacements,  and  stresses  on  the 
monitoring  surfaces  presented  in  Appendix  B  show  arrivals  at 
approximately  2.5  seconds  which  are  probably  due  to  freezing 
the  crater  motions. 
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GROUND  MOTION  RESULTS 


A  receiver  in  the  ground  close  to  the  free  surface 
tends  to  experience  two  distinct  arrivals  from  a  nearby  sur¬ 
face  burst.  These  correspond  to  the  airblast  induced  shock, 
propagating  downward  from  the  free  surface  and  the  direct  in¬ 
duced  shock  propagating  roughly  spherically  from  ground  lero. 
The  relative  magnitudes  and  order  of  arrival  of  each  depends 
upon  the  location  of  the  receiver  relative  to  both  the  free 
surface  and  ground  zero. 

Figure  4.9  shows  the  calculated  pressure  versus  time 
monitored  at  10  meters  below  initial  ground  at  a  horizontal 
range  of  80  meters.  Figure  4.10,  the  velocity  at  that  location, 
and  Figure  4.11,  the  displacement  there.  These  plots  were 
made  after  completion  of  the  calculation  by  accessing  the 
available  computer  dump  tapes  widely  spaced  in  time  compared 
to  the  time  step  of  the  calculation.  This  imolies  that  where 
gradients  are  steep,  peak  values  in  pressure  or  velocity  most- 
likely  were  not  sampled.  Nevertheless,  much  useful  in fo-ma t : on 
may  be  extracted  from  plots  of  this  type. 

The  first  peaks  in  Figures  4.9  and  4.10  at  approximately 
12  msec  are  from  the  airblast  induced  motion.  As  anticipated, 
the  velocity  is  primarily  directed  downward  into  the  ground. 

The  second  peak  at  about  30  msec  is  a  result  of  the  direct 
ground  coupling  from  the  surface  burst.  At  this  location,  it 
enhances  the  radial  (horizontal)  component  of  the  motion  and 
reverses  the  vertical  motion.  Geometry  tells  us  that  the  en¬ 
hancement  of  the  horizontal  motion  is  a  characteristic  of  this 
shock  at  all  receiver  locations  (except  directly  belov.  the 
burst  point).  The  upward  reversal  of  vertical  motion  in  local¬ 
ized  very  near  to  the  free  surface.  fit  300  msec,  Figure  4.13 
shows  a  vertical  displacement  of  approximately  5  meters  and  a 
horizontal  displacement  just  slightly  smaller.  This  station 
is  located  near  the  edge  but  within  the  calculated  crater. 
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Figure  4.9  Pressure  versus  time  at  a  horizontal  ranqe 
of  80  meters  and  a  deoth  of  10  meters. 
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Figure  4.10  Horizontal  radial  and  vertical  cowner.ts 
(Z)  of  velocity  at  a  horizontal  ranae  oc 
30  meters  and  a  depth  of  10  meters. 
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igure  4.11  Displacement  components  at  a  horizontal 
of  80  meters  and  a  depth  of  10  meters. 


Figures  4.12,  4.13,  and  4.14  show  respectively  the 
pressure,  velocity  and  displacement  components  versus  time 
at  the  same  range  of  80  meters  but  at  a  depth  of  35  meters. 
The  two  pressure  peaks  are  closer  together  in  time  due  to  the 
later  arrival  of  the  airblast  signal  from  tne  free  surface. 

At  this  depth,  Figure  4.13  indicates  the  vertical  motion  is 
downward  from  the  direct  induced  as  well.  Again,  the  primary 
motion  due  to  this  second  arrival  is  outward  horizontal.  Al¬ 
though  this  station  is  only  35  meters  from  the  free  surface, 
the  upward  vertical  displacement  and  velocity  at  300  msec  are 
approximately  7  times  smaller  than  at  a  depth  of  10  meters. 

Figures  4.15,  4.16,  and  4.17  are  similar  plots  at  the 
same  range  and  a  depth  of  75  meters.  At  this  depth,  the  two 
signal  arrivals  are  close  enough  that  it  is  impossible  to 
distinguish  one  from  the  other  on  the  pressure  or  velocity 
plots.  At  this  station  the  vertical  displacement  is  down 
into  the  ground.  Note  on  Figure  4.15  the  straight  line  from 
the  origin  to  about  40  msec.  This  is  plotting  error  which 
appears  on  many  of  these  Figures.  The  pressure  should  tc- 
zero  out  to  a  time  of  approximately  30  msec,  but  since  the 
dump  tapes  are  accessed  at  wide  intervals  o f  time  rather  than 
saving  data  at  each  calculational  cycle,  the  plot  code  has 
simply  connected  the  origin  to  the  first  accessed  data. 

Plots  of  pressure,  velocities,  and  displacements  oed 
to  300  msec  have  been  made  at  many  locations  in  order  to 
examine  carefully  the  nonlinear  ground  motion  from  the  1 -MT 
surface  burst.  Appendix  A  contains  a  representative  subset 
of  these  plots.  Included  are  plots  at  horizontal  ranees  cf 
100,  120,  160,  200  and  2  40  meters  and  at  depths  of  20,  0 , 

and  75  meters  below  the  free  surface.  At  those  depths,  i . 
plots  show  the  airblast.  motion  arrivinc  before  thr  direr- 
induced  motion.  Of  course  at;  some  greater  depth,  deter  ! :  r.  : 
on  horizontal  range,  the  direct  induced  arrival  will  bo  ;  .  -  * 
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Figure  4.15  Pressure  versus  time  at  a  ranae  oe  80  meter 
and  a  depth  of  75  meters. 
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Figure  4.16  Velocities  versus  time  at  a  range  of  80  meters 
and  a  depth  of  75  meters. 
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Figure  4.17  Displacements  versus  time  at  a  range 
meters  and  a  depth  of  75  meters. 


The  airblast  has  propagated  out  from  ground  zero  at 
speeds  initially  much  greater  than  the  elastic  wave  speed  in 
the  ground;  however,  it  rapidly  attenuated  and  slowed  down. 

At  633  meters,  the  location  of  the  CRAM-SAGE  boundary,  the 
airblast  arrived  at  approximately  0.15  sec.  with  a  peak  pres¬ 
sure  of  only  about  30  bars.  An  average  airblast  wave  speed 
from  the  burst  point  to  the  CRAM- SAGE  boundary  is  more  than 
twice  the  elastic  P-wave  speed  of  2000  meters  per  second. 
However,  the  local  airblast  wave  speed  is  already  10  percent 
less  than  P-wave  speed  at  this  location  and  rapidly  decreasing 
with  range.  This  implies  that  at  greater  ranges,  the  ground 
motion  due  to  the  airblast  at  633  meters  will  arrive  before 
the  airblast  at  the  greater  range.  Thus  the  first  motion 
near  the  free  surface  will  be  upward  out  of  the  ground, 
followed  by  a  downward  directed  motion  due  to  the  arrival  of 
the  airblast. 

Figures  4.18  and  4.19  show  the  elastic  ground  motions 
saved  at  two  monitoring  stations  8  meters  from  the  free  sur¬ 
face;  the  first  at  a  range  of  712  meters  and  the  second  at 
815  meters.  The  plots  show  both  components  of  displacements 
and  velocities  and  three  stress  components  versus  time  out 
to  3.26  sec.  Note  that  in  these  plots  the  positive  vertical 
direction  is  into  the  ground.  This  notation  is  unfortunately 
opposite  to  plots  shown  earlier  but  does  represent  the  way 
the  code  coordinates  were  oriented,  the  earlier  plots  having 
been  altered  in  order  to  show  the  free  surface  at  the  top  of 
the  grid  for  clarity  rather  than  at  the  bottom.  Thus,  as 
discussed  before,  these  plots  do  show  the  first  vertical 
motions  negative  (out  of  tne  ground).  Compressive  stress 
components  in  these  plots  are  shown  negative.  All  normal 
stresses ' shown  are  relative  to  a  hydrostatic  overburden 
pressure . 

Appendix  B  contains  a  representative  sampling  (20 
plots)  of  the  elastic  ground  motion  monitored  at  the 
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Displacements,  velocities  and  stresses  versus 
time  8  meters  from  the  free  surface  at  a  ranqi 
of  712  meters. 
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Figure  4.19 


Displacements,  velocities,  and  stresses  versus 
time  8  meters  from  the  free  surface  at  a  range 
of  815  meters. 


cylindrical  surfaces  in  SAGE  shown  in  Figure  4.6.  Since 
stresses  in  CRAM  include  the  hydrostatic  overburden  pres¬ 
sure  and  those  in  SAGE  do  not,  to  avoid  confusion,  we  have 
not  published  plots  from  CRAM  stations.  The  CRAM  results 
are  quite  similar  to  those  of  Appendix  B.  Since  all  of  the 
data  along  the  monitoring  surfaces  in  both  CRAM  and  SAGE 
are  elastic,  taking  data  from  any  of  the  three  horizontal 
surfaces  together  with  data  from  any  of  the  three  verticals 
should  give  the  same  results  when  propagated  to  the  far  field. 
Our  first  choice  for  the  cylindrical  monitoring  surface  was 
the  middle  horizontal  and  vertical  monitoring  lines  of  Figure 
4.6  (the  first  lines  in  SAGE). 

The  results  of  Appendix  B  show  motions  at  elastic 
stations  deep  below  ground  zero  to  have  virtually  ceased  by 
1.50  sec.  out  to  ranges  of  approximately  400  meters.  However, 
a  long-term  ground  roll  can  be  seen  on  the  plots  at  all 
depths  at  larger  ranges.  The  elastic  motion  due  to  airblast 
acting  on  the  free  surface  at  ranges  beyond  a  monitoring  sur¬ 
face  will  be  subtracted  from  these  data  before  proceed!:'::  with 
the  analytical  propagation  of  these  results  to  ranges  of  in¬ 
terest  which  is  the  primary  purpose  of  this  effort. 

Other  more  near  source  results  may  be  of  interest  to 
the  deep  basing,  MX  and  cratering  communities.  ’•’i.gurc  4.. 10 
is  a  plot  of  peak  pressure  versus  depth  directly  below  ground 
zero.  The  overlay  from  STREAK  to  CRAM  occurred  at  a  pressure 
level  of  approximately  40  kbar.  Figures  4.21  and  4.2?  show 
on  two  different  scales  the  contours  of  maximum  pressure  ever 
seen  at  given  locations  in  the  mesh.  Ail  contours  are  rela¬ 
tive  to  the  hydrostatic  overburden.  Since  these  contours  were 
made  by  accessing  computer  dump  tapes  which  are  spaced  many 
time  steps  apart,  they  are  quite  rough  and  probably  miss  some 
of  the  p.eaks.  Both  plots  show  the  near  surf  ice  in*  or  act  ion 
region  between  the  direct  induced  spherical  shock  and  th" 
more  planar  airblast  induced  shock.  At  depth,  the  shape  : s 
roughly  spherical. 
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Figure  4.21 


Contours  of  maximum  pressure  penetration 
showing  immediate  source  reaior. . 
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Figure  4.22. 


Contours  of  maximum  pressure  penetration 
showing  stress  levels  less  than  3  kbars. 
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Figure  4.3  gave  the  pressure  contours  in  a  portion  o; 
the  mesh  at  a  time  of  10.15  msec.  (As  dist inouished  from 
Figures  4.21  and  4.22  which  gave  peak  pressure  contours  whic.a 
have  no  times  associated  with  them.)  Figure  4.23  shows 
similar  pressure  contours  at  26.52  msec  indicating  the 
spherical  shock  propagation  and  attenuation  with  time.  Also 
shown  are  the  locations  of  interaction  with  the  planar  shock. 
In  Figure  4.23,  the  20  bar  (A)  contour  at  a  depth  of  approxi¬ 
mately  100  meters  indicates  the  hydrostatic  overburden  pres¬ 
sure  at  that  depth.  Figure  4.24  gives  the  contours  at  99.39 
msec.  Note  the  different  contour  scales  and  the  50  bar  (B) 
overburden  contour.  Both  Figures  4.23  and  4.24  show  only  a 
portion  of  the  CRAM  grid.  Figure  4.25  gives  these  contours 
over  the  entire  grid  at  647  msec,  late  in  the  calculation, 
long  after  the  shock  waves  have  passed  into  the  SAGE  grid. 

By  this  time,  the  pressures  are  beginning  to  readiust  to 
their  final  values. 

4.3  CRATER  FORMATION 

Estimates  of  crater  size  have  been  made  from  the  cal¬ 
culated  velocity  field  and  the  amount  of  tensile  crackinc  at 
late  times.  The  extent  and  orientation  of  tensile  cracks 
are  shown  in  the  CRAM  grid  plots  of  Section  4.1.  There, 

Figure  4.5  showed  tensile  cracking  for  the  first  time  at- 
51.05  msec.  These  cracks,  located  at  a  range  of  approxi¬ 
mately  100  meters  and  within  40  meters  of  the  free  surface, 
appear  to  be  oriented  toward  the  upper  left  of  the  plot. 

Figure  4.7  showed  extensive  tensile  cracking  all  around  the 
melted  blowoff  materials  at  321.91  nsec  in  addition  to  cracks 
near  the  free  surface  at  a  range  of  120  to  160  meters.  Cracks 
are  indicated  in  all  three  principal  directions .  These  cracks 
were  more  visible  in  the  rezoned  grid  a*:  64~  msec  chew n  >n 
F igure  4.8. 
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F igure  4.23. 


Pressure  contours  at  26.52  msec. 
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A  good  estimate  of  the  growth  of  the  crater  with  time 
may  be  made  by  examining  the  velocity  vectors.  Figures  4.26 
through  4.30  show  these  velocity  vectors  at  times  of  26.52, 
51.05,  99.39,  321.91  and  647  msec  respectively.  The  location 
of  the  motion  is  given  by  the  tail  of  the  vector.  Note  the 
change  in  vector  scale  from  plot  to  plot  given  at  the  top  of 
these  Figures. 

In  all  plots,  the  largest  vectors  are  in  the  cratering 
region.  (In  Figures  4.29  and  4.30,  velocities  in  the  melt 
region  have  not  been  plotted.)  Of  particular  importance  are 
the  upward  directed  velocity  vectors  near  the  free  surface 
at  the  later  times  which  indicate  the  possibility  of  those 
masses  being  ejected  from  the  crater. 

Only  those  zones  (masses)  having  sufficient  upward 
directed  kinetic  energy  to  rise  above  the  ground  surface  may 
be  removed  from  the  grid  to  form  the  crater.  Ballistic 
trajectory  calculations  determine  which  of  these  masses 
fall  back  into  the  crater  and  which  are  ejected  to  form  the 
lip  of  the  crater.  At  a  time  of  755  msec,  ballistic  calcu¬ 
lations  were  made  to  determine  the  crater  size  using  the 
kinetic  energy  escape  criteria  alone  and  also  combining  this 
criteria  with  various  tensile  cracking  criteria.  We  found 
that  requiring  the  zone  to  be  either  melted  or  cracked  for 
ejection  from  the  crater  gave  only  one  fewer  zone  ejected 
when  compared  with  the  kinetic  energy  escape  criterion. 

Based  on  this  near  agreement,  we  decided  not  to  use  any 
tensile  cracking  criteria  in  determining  the  crater  dimensions. 

Ballistic  ejecta  calculations  were  made  at  different 
times  from  0.75  sec.  out  to  approximately  1.25  sec.  Very 
little  variation  in  the  calculated  crater  volumes  (20  percent) 
were  seen.  Figure  4.31a  shows  the  CRAM  grid  at  1.17  sec  with 
the  ejected  zones  removed.  Figure  4.31b  shows  (to  the  same 
scale)  a  histogram  of  the  calculated  locations  of  the  materials 
ballistically  ejected  from  the  crater  at  this  time. 
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Figure  4.28.  Velocity  vectors  at  99.39  msec 
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Figure  4.29.  Velocity  vectors  of  321.91  msec. 
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Figure  4.30,  Velocity  vectors  at  647  msec. 
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Figure  4.31b.  Histogram  of  ejected  material. 
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Approximate  crater  dimensions  are  a  radius  of  105  meters, 
a  depth  of  90  meters,  and  a  volume  of  1.04  x  10^  m^ . 

The  calculated  crater  volume  at  1.17  sec.  corresponds 
to  a  cratering  efficiency  of  36.7  cubic  feet  per  ton  of 
explosive  energy  plus  or  minus  20  percent.  Cooper  (1976) 
has  empirically  estimated  a  cratering  efficiency  of  100  cubic 
feet  per  ton  for  wet,  soft  rock  and  200  cubic  feet  per  ton 
for  wet  soil.  Both  estimates  are  based  on  extrapolating  the 
results  of  the  Pacific  nuclear  tests  in  saturated  coral  to 
other  geologies.  The  Pacific  craters  are  unusual  in  both 
size  and  shape.  On  a  scaled  basis  they  are  considerably 
broader  and  shallower  than  the  bowl  shaped  craters  encountered 
elsewhere,  both  experimentally  and  in  calculations  of  the 
type  presented  here.  A  considerable  segment  of  the  cratering 
community  hypothesize  that  the  shape  and  size  of  the  Pacific 
craters  were  due  to  a  late  time  collapse  of  large  voids  in 
the  unique  structure  of  coral.  If  this  hypothesis  proves 
correct,  the  generic  number  for  wet  soil,  the  material  cal¬ 
culated  here,  scaled  from  the  Pacific  tests,  would  be  more 
in  line  with  the  calculations. 
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APPENDIX  A 

NONLINL,  R  GROUND  MOTION 


Plots  of  pressure,  velocities,  and  displacements 
versus  time  out  to  300  msec  at  a  series  of  locations  near  the 
free  surface  are  presented  here.  Included  are  horizontal 
ranges  (R)  of  100,  120,  160,  200  and  240  meters  and  depths  (Z) 
below  the  free  surface  of  20,  50  and  75  meters.  The  station 
location  is  given  in  the  upper  right  hand  corner  of  each  plot. 
Some  plots  show  a  straight  line  connecting  the  origin  to  the 
first  active  piece  of  data,  rather  than  a  set  of  zeroes.  This 
is  a  plotting  error  due  to  the  method  used  in  generating  these 
plots;  i.e.,  accessing  computer  dump  tapes  after  completion 
of  the  calculation.  Because  only  relatively  large  intervals 
of  time  were  sampled,  peak  values  of  pressure  and  velocity  may 
have  been  missed  in  these  plots.  For  this  same  reason,  indi¬ 
vidual  plots  have  less  wiggles  than  if  data  were  available  for 
each  computational  time  step. 
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APPENDIX  B 


ELASTIC  GROUND  MOTIONS 


Appendix  B  contains  plots  of  calculated  ground  motions 
at  20  stations  on  cylindrical  monitoring  surfaces  in  the  linear 
elastic  SAGE  code.  These  plots  are  a  representative  sampling 
of  the  elastic  ground  motion  data  which  includes  269  stations 
on  the  monitoring  surfaces  in  SAGE  as  well  as  98  stations  out¬ 
side  of  the  elastic  radius  in  CRAM.  The  data  at  each  station 
consist  of  horizontal  (X)  and  vertical  (Y)  components  of  dis¬ 
placement  and  velocity,  horizontal  and  vertical  components  of 
normal  stress  and  shear  stress.  Normal  stresses  are  relative 
to  a  hydrostatic  overburden  pressure,  a  negative  stress 
indicating  compression.  The  positive  Y  direction  on  these 
plots  is  down  into  the  ground. 
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